Late Arriving Particles in Cosmic Ray Air Showers and 
AGASA's Determination of UHECR Energies 
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We give the first detailed study of the arrival time distribution of nucleons in UHECR air showers. 
We analyze in detail the influence of late arriving particles on the energy determination of the 
AGASA experiment, as well as how the arrival time distribution changes with distance from shower 
core. Our calculations are consistent with experimental observations of the AGASA group. Crucial 
to obtaining agreement, is the correct implementation of the energy loss for low-energy protons. We 
confirm AGASA's estimation of the error in their energy determination associated with late-arriving 
particles, assuming primary protons. 
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I. INTRODUCTION 



The results of the AGASA collaboration have been fol- 
lowed with great interest throughout the years because 
they had the largest exposure for highest energy cosmic 
rays. As of now, AGASA reports about 10 cosmic ray 
events with an energy greater than 10 20 eVpJ, well be- 
yond the predicted GZK cutoff due to interaction of par- 
ticles with the cosmic microwave background, assuming 
uniformly distributed sources. By contrast, the HiRes 
collaboration reports that their monocular results are 
consistent with a GZK cutoff 2] . The statistics are not 
yet sufficient to decide whether this is a significant dis- 
crepancy or not, but the question of whether there is a 
systematic shift between the relative energy normaliza- 
tions of AGASA and HiRes has assumed special impor- 
tance. AGASA is an extensive air shower array detector 
while HiRes uses the air fluorescence technique. 

In this paper we investigate the contribution of nucle- 
ons to AGASA's energy determination, which has been 
neglected in previous analyses. The energy determina- 
tion in AGASA is based on the particle density, measured 
with plastic scintillator detectors. The reconstructed 
value of the signal at a distance of 600m from the core 
is used as an energy estimator, since this is considered 
to be rather insensitive to fluctuations and to not de- 
pend too much on the air shower model. The density of 
nucleons at 600m is not negligible. Neutrons (unlike pro- 
tons) have a small energy deposit in plastic scintillator 
but they scatter elastically with protons in the scintil- 
lator. The arrival-time-delay distribution of neutrons is 
very broad and, due to AGASA's technique for recording 
the signal, late arriving particles have an exponentially 
increased weight as is explained below. This is the moti- 
vation for our scrutiny of the issue of late arriving parti- 
cles in general and nucleons in particular. HiRes' energy 
measurement entails an entirely different set of issues and 
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FIG. 1: Lateral distribution of particle density as a function of 
distance from the core in a 1 x 10 19 eV proton-induced shower. 



systematic uncertainties which we do not address in this 
paper. 



II. GROUND NUCLEONS 

The air-showers analyzed in this paper are simulated 
with the SENECA model, as introduced in 0. The ba- 
sic features of this model are the high energy hadronic 
model QGSJET01 0, the electro- magnetic shower model 
EGS4 5] , and different choices of the low-energy hadronic 
model: GHEISHA 6], which was the default option in 
the CORSIKA model, G-FLUKA 7] and GCALOR @ 
used in the framework of the GEANT 3.210, pack- 
age for detector simulation. The model also has some 
new simulation techniques which speed up the computa- 
tion of air showers considerably, by using the approach of 
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FIG. 2: Arrival time distribution for particles at 600 m from FIG. 3: How late arriving particles depositing 10% of the 
the shower-core. Times are with respect to the arrival of the initial energy deposit can fake a signal 1.7 times higher, 
shower- core. 



cascade equations as described in detail in ref. [lfj • The 
pure physics content (hadronic and electromagnetic mod- 
elling) using GHEISHA as low-ener gy h adronic model, is 
identical to the one of CORSIKA [11|, and the results 
are in good agreement with this model. 

Figure^shows the averaged lateral distribution of par- 
ticle densities as a function of the distance from the 
shower core, generated by 1 x 10 19 eV vertical pro- 
ton showers at an altitude 667m above sea level, using 
GCALOR as the low-energy hadronic model. All par- 
ticles have been plotted which have kinetic energies l 
1 MeV. Neutrons become quite prominent at large dis- 
tances from the shower-core, unlike protons which are 
more readily absorbed in the atmosphere due to ioniza- 
tion energy loss. 

Electromagnetic particles and muons arrive within 
about a micro-second, whereas nucleons can be retarded 
by many tens of microseconds due to multiple scatter- 
ings. The arrival time distribution of particles of various 
species, 600 m from the core, is shown in Fig. |3 Ar- 
rival times are measured with respect to the arrival of 
the shower core. 



III. PARTICLE MEASUREMENT IN AGASA 

A particle of a specified type, kinetic energy and inci- 
dence angle, gives a certain energy deposit in AGASA's 
5 cm thick plastic scintillator. When a charged particle 
passes through the dectector, the signal from the photons 
generated in the scintillator and detected by the photo- 
multiplier is processed to produce a signal which decays 
exponentially, with a time constant of 10/is. The time du- 
ration in which the signal has an amplitude greater than 
a given discrimination level is called the pulse-witdh. 
Recording the pulse-width rather than the signal itself 
has the advantage that a large dynamic range of energy 



deposit is measurable, since the pulse-width is propor- 
tional to the logarithm of the energy deposit. 

An important assumption is that all particles arrive in 
a time much shorter than the decay constant of the signal. 
The bulk of particles arrive in less than a micro-second, 
but neutrons can scatter for a much longer time since 
they do not lose energy by ionization of air molecules, as 
is reflected in Fig. 2 The effect of a particle arriving 20/ts 
late, with an energy deposit 10% as large as the original 
signal, is shown in Fig. 01 Even though the late energy 
deposit in this particular example is small, the time delay 
increases its importance by a factor of exp(20/xs/10/xs) = 
7.4 and the inferred total energy deposit is overestimated 
by roughly a factor of 1.7 . 

One can easily calculate the effect of this time delay 
for the measurement of particle densities. It increases the 
signal by the factor 

ft: dJ ^eM(t-to)/10ys)dt 

J to dt UL 

where to is the time when the detector triggered and t\ 
is the time when the signal falls below the threshold or 
128/xs, whichever is less, t^ represents the arrival time of 
the last particle, not necessarily equal to t\. Thus / is the 
factor by which the measured signal has to be reduced in 
order to obtain the true energy deposit in a detector. 

Figure 0] shows / as a function of radius for the total 
signal of electrons, muons and photons, using the shower 
simulation described in section^] The factor / increases 
as a function of the radius due to greater spreading of the 
arrival times at large distance. At very large distances / 
decreases again, when the particle number in a detector 
falls to of order unity. If a single particle arrives at a 
detector, / = 1 as is evident from replacing dE^cp/dt in 
eq. (JIJ with a delta function at to- However / can also 
be less than 1, since in general the arrival time ti can 
be larger than the time t\ when the pulse falls below the 
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threshold. If for example a second particles arrives af- 
ter t\ the AGASA acquisition system generates a second 
signal but it is ignored in the analysis. Effects arising 
when the particle number is of order one are a function 
of the area of the detector and primary energy as well 
as distance to the core. The AGASA collaboration mea- 
sured the effect of the arrival time distribution, in a 30 
m 2 detector with time-resolution as described in ref. 0. 



IV. EFFECT OF NEUTRONS AND PROTONS 
IN A PLASTIC SCINTILLATOR 

The energy-deposit of neutrons in a 5cm-plastic scintil- 
lator is calculated with the GEANT-package |j and the 
Gcalor interface @ relevant for low energy neutrons be- 
low 10 GeV. Above 10 GeV, routines from the FLUKA 
package are called automatically Q- The geometrical 
setup consists of a 1.5m x 1.5m x 5cm plastic scintillator 
in a 1.6mm thick iron box which is placed in a farmer's 
storage house made with 0.4mm thick iron plates. We 
have checked that we can reproduce the energy deposit 
spectra used by the AGASA-collaboration for electrons, 
positrons, muons and photons |12| . Protons deposit en- 
ergy in a plastic scintillator since they are charged par- 
ticles. The main process is ionization. Neutrons affect 
the scintillator by scattering elastically with the protons 
of hydrogen abundantly available in plastic scintillators; 
the recoil protons then deposit energy |l3j. 

The energy deposit spectra of neutrons and protons are 
shown in figure |S] for particles arriving at different angles. 
One sees the increasing energy deposit due to the longer 
path length in the material scaling with 1 / cos(#) . 

If the energy deposit in a scintillator is locally concen- 
trated, a saturation of the scintillation response is ob- 
served. The impact of this on the scintillation signal 
generated by neutrons and protons is described by Birks' 
law Q, and is shown in Fig. Throughout our calcu- 
lations, the effect of Birks' law has been included, using 
the parameters proposed in the GEANT manual which 
are consistent with values proposed in |14j . The result is 
expressed by quoting an effective energy deposit, which 
would give the equivalent scintillator response in the ab- 
sence of this effect. 



V. CALIBRATION OF THE THEORY 

AGASA calibrates their detectors with the ambient 
signal from atmospheric muons, also called "omnidirec- 
tional" muons, which are the dominant source of signals 
in individual detectors. Because these muons are close 
to being minimum-ionizing, the energy deposited by a 
single muon depends primarily on its angle (tracklength) 
and depends only weakly on its energy. Therefore the 
pulsewidth distribution of an individual detector is well 
understood and suitable for calibrating the scintillator 
response. We calibrate the theory in the same way. 
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FIG. 4: The effect of time delay on the measurement of the 
total energy deposit of electrons, muons and photons. 
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FIG. 5: Energy deposit in a 5cm plastic scintillator of verti- 
cally arriving neutrons (top) and protons (bottom) of different 
energies. 
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FIG. 6: The effect of Birks' law on the effective energy depo- 
sition spectrum of neutrons. 
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FIG. 7: The fitted distributions of atmospheric muons. 

iFrom data of refs. and Caprice 94 data [t7| . 

a functional fit of the muon spectrum dD{ja,6) / dp is 
obtained; data and fit are shown in Fig. [7| Then a 
Monte Carlo simulation is performed, in which muons are 
produced with this distribution and propagated through 
our detector simulation. Here, it is crucial to apply 
a smearing function to the energy deposit to repro- 
duce the experimental shape of the pulse width distri- 
bution. We have found that a suitable smearing func- 
tion is P{x) ~ exp(— 0.5 la 2 (xb)/a 2 )/x, which can be 
obtained by taking the exponential of a random num- 
ber sampled from a Gaussian distribution around zero 
and with width a. Empirically, this function imitates 
Landau fluctuations in the photo-multiplier tube. The 
parameter b is adjusted such that the mean is equal to 
one: J xP(x)dx/ J P(x)dx — 1. The values a = 0.3 and 
b = 1.045 proved to describe best the shape of the PWD 
observed by AGASA. The quality of this fit is shown in 



FIG. 8: Theory calibration according to the experiment. 
Each channel corresponds to 500ns. AGASA calibrates with 
the peak value ii whereas we calibrate the threshold of the 
signal (energy deposit) such that the peak is at t\ = 10/is. 

Fig.0 

The next step in calibration is to set the absolute nor- 
malization of the signal. To convert AGASA's calibra- 
tion to an energy scale we need to return to a more de- 
tailed discussion of their data recording procedure, which 
was described schematically in section ILTT1 If the signal 
threshold is set at Vthr , then the time-above-threshold or 
pulse- width, tp, is implicitly given by [jj] 

V thr =VeM~tpM , (2) 

where V is the pulse-height of the signal, proportional to 
the scintillation response and therefore to energy deposit. 
A priori, V and Vthr are measured in an arbitrary unit 
relating to the gain of the PMTs. 

By defining t\ as the pulse-width of a single particle 
with pulse-height V\ one obtains 

V thl = Vi expt-tx/r) , (3) 

and 

^=^ = exp((tp-ti)/r) (4) 

as the number of particles for a arbitrary pulse- width ip . 
One can consider as some effective number of particles, 
since it comprises different particle types. AGASA sets 
their gain and threshold so that the peak of the pulse- 
width distribution of atmospheric muons occurs roughly 
in channel 20, as shown in Fig. |S| This corresponds to 
t\ fa t = 10fj,s because one channel is 500ns. The precise 
calibration is then given by the peak value t\. 

In our case the signal is just the energy deposit E in 
the plastic scintillator, so V — > E applies to all formulas 
above. We calibrate the threshold E'thr such that the 
peak of the pulse-width distribution of omnidirectional 
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FIG. 9: LDF for a 10 19 eV vertical proton induced shower 
using QGS JET01 /GHEISHA. 



muons is exactly at ti = r = 10/is. According to J2J, 
the energy deposit of an effective particle is then E\ = 
£"thr/ exp(— 1) and for arbitrary tp the effective number 
of particles is just 



N = exp(i P /r - 1) 



(5) 



Our simulation shows that the peak in the energy de- 
position of the omnidirectional muons corresponds to 
Ei = 0.011 GeV per particle. 



VI. LATERAL DISTRIBUTION FUNCTIONS 

The AGASA collaboration finds experimentally that 
the lateral distribution of particle density (as defined via 
equation can be described by the empirical function: 



R 

Rm 



-(rj-a) 



1 



(6) 

where the parameters a = 1.2 and 6 — 0.6 are fixed, 
Rm — 91.6m is the Molicre radius two radiation lengths 
above the altitude of the Akeno observatory, and the 
slope 77 at distances large compared to Rm is an experi- 
mentally determined function of the shower zenith angle 
6: r) = 3.84- 2.15(sec<9 - 1) Q ■ 

By fitting "data" from their Monte Carlo shower and 
detector simulation to the lateral dependence of equa- 
tion ©, AGASA obtained the conversion formula (for- 
mula (1) from Ref.|lJ ) 



E = 2.03 x 10 17 S'(600m)eV, 



(7) 



using calculations for the Akeno observatory at the alti- 
tude 900 m. It was realized later that the average altitude 
of the whole AGASA array is actually 667 m, significantly 
different from that of the Akeno sector only. However 
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FIG. 10: LDF scaled by the AGASA empirical function for 
different low energy models and with the effects of the shower 
front and late nucleons. 



when the AGASA analysis is corrected for the change 
in average altitude and still other effects (mainly shower 
front thickness and delayed particles - see table 2 from 
) , AGASA finds a final conversion formula which is co- 
incidentally the same as their original one, J7J|. Below, 
when we refer to the AGASA LDF, we mean equation 
(6) normalized by 0. 
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We now compare the results of our calculations for 
various cases, to AGASA's empirical results. We use an 
average altitude 667 m and the procedure of Sec. to 
obtain the signal as a function of energy deposit. Fig. [5] 
shows our predicted LDF for an average 10 19 eV proton- 
induced vertical shower compared to the the AGASA 
empirical function @. To quantify the importance of 
including the shower-front thickness and the contribu- 
tion of nucleons, the figure shows three cases, all us- 
ing QGS JET01 / GHEISHA as the hadronic models. The 
curves are labeled E dep , LDF and LDF+nucleons, and 
are computed as follows: 

• Edep- Only the energy deposit of elec- 
trons/positrons, muons, and photons is included; 
particles are assumed to arrive simultaneously. 

• LDF: Same particles contributing as above, but 
with non-trivial distribution of shower front thick- 
ness taken into account. 

• LDF+nucleons: the above, and additionally the ef- 
fect of nucleons and their time delay taken into ac- 
count. 

One observes that the combination of hadronic mod- 
els QGS JET01 /GHEISHA, gives a slightly too flat slope 
of the LDF in all cases. This has already been found 
in ref. [l^ based on "Edep" alone. We see here that 
including the shower front thickness and the nucleons in- 
creases the discrepancy. This can be seen better in Fig. 
HUIwhere the results are shown as a ratio to the AGASA 
LDF. The lower panels of this figure show the LDF as 
obtained by choosing G-FLUKA and GCALOR as low 
energy hadronic model instead of GHEISHA. One sees 
that these two give a better description of the AGASA 
LDF, though still dropping slightly more slowly with dis- 
tance. That the tails of LDFs can be sensitive to the low 
energy hadronic model has already been shown in ref. 

The values of the scaled LDFs at 600 m in Fig. QJJ al- 
low us to infer the deviation of AGASA's energy estima- 
tion equation (JJJ from that implied by our simulations. 
Table [|] shows the ratio of our predicted signal to that 
"predicted" by AGASA using J7J) as conversion factor, 
with QGSJET01 for the high energy model in each case. 
Thus for a given observed signal, if GHEISHA were the 
correct low-energy hadronic model AGASA would overes- 
timate the energy of a vertical proton by about 5%, while 
if G-FLUKA or GCALOR correctly describe the low en- 
ergy interactions, the AGASA energy estimate would be 
about 2% too low, when all the effects of late-arrivers 
and shower-front thickness are included. 



VII. COMPARISON OF PREDICTED AND 
OBSERVED ARRIVAL TIME DISTRIBUTIONS 

As discussed in Sec. |Hj given AGASA's method of 
recording the signal, late arriving particles such as nu- 



Model deviation from £J 

QGS JET01 /GHEISHA 1.049 
QGS JET01/ G-FLUKA 0.983 
"OGSJETOl/GCALOR I 0.978 

TABLE I: The deviation of the energy estimation obtained 
in these simulations. 
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FIG. 11: Over-estimation factor due to late arriving nucleons 
as a function of radius, for 10 19 eV proton-induced showers. 



cleons could in principle distort the determination of the 
primary energy; Fig. 1111 shows the simulated value of 
the overestimation factor / due to late arriving nucle- 
ons as a function of distance to the core. However as we 
saw in the previous section, the impact of late arrivers is 
compensated by other effects which had originally been 
neglected in AGASA's calculations, so that we find the 
average AGASA energy determination to be quite good, 
assuming the primaries are protons and QGSJET01 is an 
adequate high energy model. The overestimation factor 
at 600m is not more than 5% for a vertical 10 19 eV pro- 
tons and we have confirmed it is the same for 5 x 10 19 eV 
primaries. 

The arrival time distribution itself has not up to now 
been critically compared to observations, which we do 
in this section. In addition to the standard detectors, 
the AGASA collaboration^] used a 30 m 2 detector in 
combination with the rest of the AGASA array to record 
the actual signal in that detector as a function of time, 
within ± 30/xs of the first particles of the shower. Ref. pj 
remarked that including nucleons in model calculations 
gives an LDF which is much too flat compared to what 
is observed experimentally. We have found that this be- 
havior can be produced by an imprecise implementation 
of energy loss. 

The Bcthe Bloch formula for the energy loss of charged 
particles diverges for values of the relativistic gamma- 
factor approaching 7=1. When particles are propagated 
over a large distance in a single step, this divergence in 
dE/dx is not properly taken into account if the Simula- 
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other than nucleons. Energy loss is accurately treated in 
electromagnetic shower codes, and when other charged 
particles such as muons and pions in the shower become 
non-relativistic they decay too quickly for their subse- 
quent energy loss to matter. All results shown in this 
paper have been computed with high precision, for all 
charged particles in the shower. 

Now we compare the observed delay time to our sim- 
ulations. The lower panel of Fig. ^| shows the fraction 
of delayed particles as a function of the distance from 
the shower-axis, calculated with high and low precision. 
The data is from the AGASA collaboration measurement 
discussed abovep] and includes isotropically distributed 
showers up to 45° incident angle. The simulations are 
done in the same way, assuming purely proton primaries. 
The fraction of delayed particles is defined as the energy 
deposit of particles arriving later than 3 /is after the first 
particles in the detector, divided by the energy deposit of 
all particles. The late arrivers in the high precision cal- 
culation are reduced by a factor of two, and the results 
are in good agreement with the data. 



VIII. SUMMARY AND CONCLUSIONS 
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FIG. 12: LDF (upper panel) and fraction of delayed particles 
(lower panel) with high precision (hp) and low precision (lp) 
calculations of energy loss of charged particles. 



tion calculates dE/dx at the starting point only and takes 
it to be the same over the total step length. If one takes 
for the step length the distance to the next interaction or 
decay, it is for non-relativistic particles often much too 
long for constant dE/dx to be a good approximation. 
For this reason the energy loss of protons, and neutrons 
coming from them as secondaries in the shower, is under- 
estimated. A simple cure for this problem is to limit the 
step length to a suitably low value Z ma x- We found that 
for l max < 10 m the results are stable, i.e., the low en- 
ergy spectra of hadrons do not change any more as l max 
is reduced further. 

The influence of the higher precision achieved can be 
seen in the upper panel of Fig. ^] where the LDF, in- 
cluding nucleons, has been calculated with high- and low- 
precision energy loss. "Low-precision" means taking the 
step length as proposed by the total cross-section of the 
interaction, usually several hundred meters at ground 
level, as done in standard simulations. It can be seen 
that the low precision energy loss calculation results in a 
significant overestimation of the LDF at large distances. 

The change in precision is not important for particles 



We have analyzed the effect of late arriving particles on 
the energy determination of the AGASA experiment. We 
find that the net effect of the arrival time spread in the 
electromagnetic and muonic components of the shower, 
and the inclusion of nucleons in the simulation, has only 
a small effect on the AGASA energy determination. It 
shifts the true energy compared to the energy as deter- 
mined without these effects, upward by a few percent 
in the models which give the best agreement with the 
shape of the lateral distribution function. This agrees 
with the correction applied by AGASA Q. The LDF 
measured by AGASA is best described with G-FLUKA 
or GCALOR as the low energy hadronic model; the cor- 
responding LDF using the GHEISHA code is too flat. 

We have demonstrated for the first time that shower 
simulations can describe the distribution of arrival times 
of late-arriving particles, primarily nucleons, as observed 
by the AGASA collaboration|lJ. When analyzing neu- 
trons and protons in air shower simulations, the energy 
loss of charged particles has to be done with care. It 
is straightforward but crucial to compute these energy 
losses accurately, otherwise the contribution of these par- 
ticles to the signal can be greatly overestimated. The 
Pierre Auger Observatory records detailed arrival time 
information in each tank, so accurately modeling the ob- 
served arrival time distributions should be a powerful 
new tool for validating shower simulations. An inter- 
esting question, left to the future, is to what extent the 
arrival time distribution can be useful for studying the 
composition of UHE cosmic rays. 
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